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Abstract 

The theory of multidimensional persistent homology was initially devel- 
oped in the discrete setting, and involved the study of simplicial complexes 
filtered through an ordering of the simplices. Later, stability properties of 
multidimensional persistence have been proved to hold when topological 
spaces are filtered by continuous functions, i.e. for continuous data. This 
paper aims to provide a bridge between the continuous setting, where sta- 
bility properties hold, and the discrete setting, where actual computations 
are carried out. More precisely, a stability preserving method is developed 
to compare rank invariants of vector functions obtained from discrete data. 
These advances confirm that multidimensional persistent homology is an ap- 
propriate tool for shape comparison in computer vision and computer graph- 
ics applications. The results are supported by numerical tests. 
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1 Introduction 

In this paper we present a discrete counterpart of the theory of persistent homol- 
ogy of vector functions that still guarantees stability properties as the continuous 
framework. The theory of multidimensional persistence was developed in the dis- 
crete setting in |[8l, and involved the study of simplicial complexes filtered through 
an ordering of the simplices. On the other hand stability properties of multidi- 
mensional persistence are proved to hold when triangulable spaces are filtered by 
continuous functions, i.e. for continuous data [|20l HI. This paper aims to be a 
bridge between the continuous setting, where stability properties hold, and the 
discrete setting, where actual computations are carried out. More precisely, we 
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develop a method to compare persistent homologies of vector functions obtained 
from discrete data. We show that in the passage from the continuous to the discrete 
framework stability is preserved. These advances support the appropriateness of 
multidimensional persistent homology for shape comparison by functions. 

The problem of comparing shapes is well-studied in computer vision and com- 
puter graphics and many algorithms have been developed for this purpose. A 
widely used scheme is to associate a shape with a shape descriptor, or a signa- 
ture, and comparing shapes by measuring dissimilarity between descriptors. An 
important class of shape descriptors, which may be called shape-from-functions 
methods, is based on the common idea of performing a topological exploration 
of the shape according to some quantitative geometric properties provided by a 
(measuring) function defined on the shape and chosen to extract shape features 
0. 

The simplest topological attribute of a space is the number of its connected 
components. A well-known mathematical tool to count the number of connected 
components is the homology group Hq. More complex topological features are 
revealed by higher homology groups. 

Persistent homology is a shape-from-functions method for shape description 
involving homology groups of any degree. The idea is to filter a space by the sub- 
level sets of the function and to analyze the homological changes of the sublevel 
sets across this filtration, due to the appearance or disappearance of topological 
attributes, such as connected components. Features with a short persistence along 
the filtration can be regarded as negligible information due to noise or very fine 
details. For application purposes, it is often sufficient to disregard the group struc- 
ture of persistent homology and retain only the rank information. This gives rise 
to the notions of rank invariant (SI, persistent Betti numbers |[T9l , size functions 

The topic has been widely studied in the case of filtrations induced by scalar 
continuous functions (i.e. one-dimensional persistence), especially in connection 
with the stability problem ffT3l[T2l [T4l[T5]| . 

This theory has been generalized to a multidimensional situation in which a 
vector- valued function characterizes the data as suggested in [fTTlfTSl . Results in 
this area are given in [[3l [H |5l |9l . This generalization is quite natural in view of 
the analogous generalization of Morse Theory ll24l . Moreover, it is motivated by 
applications where data are more completely described by more than one function 
(e.g., curvature and torsion for space curves). 

The passage from scalar to vector-valued functions presents new challenges. 
To begin with, critical points are no longer isolated even in non-degenerate situ- 
ations [[TtI . Although the relevant points for persistent homology of vector func- 
tions are a subset of the critical points, precisely the Pareto critical points, these 
are still non-isolated [jUJ. For example, in the case of the sphere + y"^ + = 1 
with the function / = (y, z), the Pareto critical points are those in the set a; = 0, 
y"^ + z"^ = 1, yz > 0. 

Another delicate issue is passing from the comparison of continuous models 
to that of discrete models. This is an essential passage, and the core of this pa- 
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per. Indeed, for two given real-world objects X and Y, modeled as triangulable 
topological spaces (e.g., manifolds), we usually only know simplicial descrip- 
tions /C and C of them, affected by approximation errors. For example, acquiring 
3D models of real-world objects for computer graphics applications needs to ac- 
count for errors due to sensor resolution, noise in the measurements, inaccuracy 
of sensor calibration fl]. Moreover, different techniques for reconstructing the 
geometry and topology of the scanned object yield different polyhedral approxi- 
mations. Analogous considerations hold for any continuous measuring functions 
/ : X — )■ M'^ and g : Y M'^, because we could only consider approximations 
: K ^ MJ' , ijj : L ^ M.'' of f and g defined on finite polyhedra. Depending on 
the context of a specific application, these functions may or may not be given by 
explicit formulas. In either case, it is legitimate to assume that we are able to com- 
pute their values on vertices of /C and C. Hence, we only know the discrete maps 
: V(/C) —7- MJ', if) : V(£) —7- M'^ which are the restrictions of (p and ip to vertices. 
Therefore a natural question is whether shape comparison by persistent homology 
of vector functions is numerically stable, i.e. whether the computation of a dis- 
tance between rank invariants of discrete models gives a good approximation of 
the ideal distance between rank invariants of continuous models. 

Our main result. Theorem 14. 5 [ gives an affirmative answer to this question. It 
states that, in the passage from continuous to discrete data, the distance between 
rank invariants does not increase, provided that stability holds for the continuous 
model. We underline that at least one stable distance between rank invariants of 
continuous vector functions exists as proved in [l9l- In order to profit from the 
stability theory in the continuous case, we give a new construction of axis-wise 
linear interpolation Lp" which is generic in the sense that its persistent homology 
is exactly equal to that of the map defined on vertices Lp. 

The paper is organized as follows. In Section |2] the necessary background 
notions concerning persistence are reviewed and put in the context of our aims. 
Section [3] starts with the description of the simplicial framework and with Exam- 
ple [3T| which is a simplicial analogue of the sphere example pointed above. The 
same example shows that, in the vector case, the linear extension of a map defined 
on vertices does not satisfy the genericity property described above. Topological 
artifacts of an interpolation method have been observed before. This phenomenon 
can be referred to as topological aliasing. Our example motivates the construction 
of our axis-wise linear interpolation. We next prove Theorem 13. 3 1 on the deforma- 
tion retraction of continuous sublevel sets of ip^ onto the simplicial sublevel sets 
of ip. We introduce the notion of homological critical value for vector functions. 
As in the sphere example, the set of critical values need not be discrete, but we 
prove in Theorem 13. 5 1 that in the case of ^p^ it has to be contained in a finite union 
of hyperplanes, thus it is a nowhere dense set and its A;-dimensional Lebesgue 
measure is zero (Corollary 13.61 ). 

Section m starts with Lemma |4n that provides an approximation of a distance 
between the rank invariants of continuous functions by that of the rank invariants 
of the corresponding axis-wise linear approximations. The genericity of y?^ allows 
us to introduce the rank invariant for ip. Although this rank invariant is defined 



3 



for a discrete function (/? and computed using only simplicial sublevel sets, it takes 
pairs of real vectors as variables, as it is in the case of the rank invariant for 
continuous functions. We show that this new rank invariant for ip is equal to that 
of If" . This allows us to derive the main result of the paper (Theorem 14 .5 1) . 

Section [5] describes an algorithm which computes an approximate matching 
distance. Our algorithm is a modification of the algorithm described in [2], adapted 
to the rank invariants. The correctness of the algorithm is guaranteed by the re- 
sults of Section m We next present tests of the algorithm performed on simplicial 
models in the case k = 2. Our tests revealed the same discrepancy as observed in 
Example |3.1[ thus providing numerical confirmation of topological aliasing. 

2 Basic notions and working assumptions 

Let us consider a triangulable topological space X (i.e., a space homeomorphic 
to the carrier of a finite simplicial complex). A filtration of X is a family = 
{Xa}aeR'' of subsets of X that are nested with respect to inclusions, that is: Xa C 
Xj3, for every a ^ /3, where a ^ 13 if and only if aj < /3j for all j = 1, 2, . . . , /c. 

Persistence is based on analyzing the homological changes occurring along 
the filtration as a varies. This analysis is carried out by considering, for a ^ 13, 
the homomorphism 

: H,{Xa) ^ H^Xp). 

induced by the inclusion map i^"-'') : Xa X^. We work with Cech homology 
with coefficients in a given field F. When each Xa, a E M'^, is triangulable, 
it reduces to simplicial homology. For simplicity of notation we write H^,{Xa) 
for the graded homology space H^{Xa] F) = {Hq{Xa] F)},^^- The choice of a 
field is only made in experimentations, the most convenient in computations being 
F = Zp, with p a prime number. Thus, for any q E "L, Hg{Xa) is a vector space 
of dimension equal to the g'th Betti number of Xa- 

The image of the map Hg(i^°''f^^) is a vector space known as the q'th persistent 
homology group of the filtration at {a, (3). It contains the homology classes of 
order q bom not later than a and still alive at (3. The dimension of this vector 
space is called a q 'th persistent Betti number. 

A rank invariant is a function that encodes the changes in the persistent Betti 
numbers as a and (3 vary. Setting 

:= {{a, (3) e M'^ X M'^ I a ^ /?}, 

where a -< /9 if and only if aj < (3j for all j = 1, 2, . . . , fc, the q'th rank invariant 
of the filtration F is the function pj: : — )■ N U {oo} defined on each pair 
{a, 13) G as the rank of the map Hq{i^°''^^). In other words, p'jr{a,(3) = 
dimimi/g(z("'^)). 

In this paper, we will use the notation pjr to refer to rank invariants of arbitrary 
order. Ultimately, the shapes of two triangulable spaces X and Y, filtered by 
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and Q, respectively, can be compared by using an (extended) distance D between 
their rank invariants pjr and pg. 

The framework described so far for general filtrations can be specialized in 
various directions. We now review the two most relevant ones for our paper. 

2.1 Persistence of sublevel set filtrations 

Given a continuous function / : X — > M''", it induces on X the so-called sublevel 
set filtration, defined as follows: 



We will call the function / a measuring fi^nction and denote the rank invariant 
associated with this filtration by pf. 

Since X is assumed to be triangulable and / is continuous, p/(«, (3) < +00 
for every a ^ /3 e M'^ (cf. 

Among all the (extended) distances D between rank invariants of filtrations, 
we confine our study to those ones that, when applied to sublevel set filtrations, 
satisfy the following stability property: 

(S) For every f,f':X^ M.'^ continuous functions, D {pj, pf) < ||/ — /'||oo- 

In ^ it has been shown that there is at least one distance between rank in- 
variants, the matching distance, that has the stability property (S). An analogous 
stability property for a distance defined between modules is presented in [22|. 

The matching distance will be used for computations in the experiments de- 
scribed in Section |5l Until then, we will not need to specify which distance D we 
are using, provided it satisfies (S). 

2.2 Persistence of simplicial complex filtrations 

We consider a simplicial complex JC consisting of closed geometric simplices and 
its carrier defined by 



The set of all vertices of /C is denoted by V(/C) or by V, if /C is clear from the 
context. For cr, r G /C, the relation r is a face of a is denoted by r < cr. For 
proper faces, we write r < a. 

In this discrete setting, we take a family {/CajagRfc of simplicial subcomplexes 
of /C, such that JCa is a subcomplex of /C/?, for every a ^ /3. As a consequence, 
their carriers are nested with respect to inclusions, that is: Ka C Kjs, yielding a 
filtration of K. 

In the next section we address the following problem: is any simplicial com- 
plex filtration induced by a suitable continuous function? 

A positive answer to this question will allow us later to transfer the stability 
property of D from the continuous to the discrete setting. 



X^ = {xeX\ f{x) ^ a}. 




(1) 
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3 From continuous to discrete vector functions 



We let ip : V(/C) — j- M.^ be a vector- valued function defined on vertices. We 
suppose that </? is a discretization of some continuous function (p : K ^ W'. 
Reciprocally, is an interpolation of ip. In this section we will simply assume 
that (p is equal to ip on vertices of K, but, of course, when it comes to comput- 
ing, one has to set bounds for the rounding error. Although in some practical 
applications of persistent homology to the analysis of discrete multidimensional 
data ip : K ^ may be explicitly known, in some other cases we do not even 
have an explicit formula for ip: we assume that such a function exists, that we 
can estimate its modulus of uniform continuity (for the sake of simplicity, say, 
its Lipschitz constant), and that we can compute the values of (p at grid points of 
arbitrary fine finite grids. 

In a discrete model, we are interested in simplicial sublevel complexes 

/Ca := {cr G /C I (p{v) :< a for all vertices v < a}. 

In Section |51 we compute the rank invariants for the discrete vector-valued 
function ip and we use this information for computing the distance between rank 
invariants for their continuous interpolations. In order to do this, we need to know 
that there exists a continuous function which is a generic interpolation of ip, in the 
sense that its rank invariant is exactly equal to that of (/?. 

In the case k = 1, that is, when ip has values in M, it can be shown that such an 
interpolation can be obtained by extending ip to each simplex cr G /C by linearity. 
We shall denote this interpolation by ip. In that case, one can show that is 
a deformation retract of K^<a, so the inclusion of one set into another induces 
an isomorphism in homology. This result belongs to "mathematical folklore": it 
is often implicitly used in computations without being proved. The arguments 
for that case are outlined in [23, Section 2.5] and Theorem 13.31 we prove in this 
section contains this result as a special case. Unfortunately, if A; > 1, this result is 
no longer true as the following example shows: 

Example 3.1 Let K be the boundary of the tetrahedron shown in FigurefH home- 
omorphic to the 2D sphere. The corresponding simplicial complex /C is made of 
all proper faces of the 3D simplex [vq, Vi, V2, v^] in M^, with vertices vq = (0, 0, 0), 
Vi = (1, 0, 0), V2 = (0, 1, 0), = (1/2, 0, 1). Hence K = \JC\ is homeomorphic 
to a 2D sphere. Let (p : K ^ M? he the restriction of the linear function Tp given 
by ^(x, y, z) = (x, z) to the four vertices. Let a G be any value chosen so 
that 1/2 < «! < 1 and a2 = 2 — 2ol\. It is easy to see that = [vq, V2]. Its ho- 
mology is trivial. Note that the set Kj^^a contains one point x on the edge [vi, v^], 
namely x = {ai, 0, 02), which closes a non-contractible path in K^^a- We have 

The discrepancy between the discrete and linear interpolated models seen in 
Example [3T| has been observed in applications to computer graphics and imaging, 
and has been recently referred to as topological aliasing. 
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Figure 1 : The tetrahedron boundary and sketches of two sublevel sets of the linear 
interpolation ip discussed in Example 13. 1[ The values taken at the displayed edge 
are critical. 



Several interesting conclusions can be derived from this. First, is not a de- 
formation retract of K^^a- This remains true if we slightly increase the value of a. 
Secondly, if we slightly decrease a, the set Ka does not change but the set Ki^^a 
becomes contractible. Hence, in the sense of Definition 13.41 presented further in 
this section, any value assumed at a point of the edge [^1,^3] is a homological 
critical value. In particular, the set of such values may be uncountable. This is 
in contrast with the one-dimensional case, where a piecewise-linear function on a 
simplicial complex must have a discrete set of critical values. 

We shall now construct a continuous function Lp^ : K M.'^ called axis-wise 
linear interpolation of p which will correct the problem encountered with the 
linear interpolation Tp in the multidimensional case. First, given any a E IC, let 

/i(cr) G M*-' be defined by 

yUj((T) = max{(pj{v) \ V isa vertex of a},j = 1, 2, . . . , /c. (2) 

Note that if r < ci, then fi^r) ^ /u(cr). 

We will use induction on the dimension m of a to define (/? : A' — )■ M''" on a 
and a point E a with the following properties: 

(a) For all x E a, (p^{x) ^ p'iwa) = l^io') ; 

(b) pP is linear on any line segment \wa, y] with y on the boundary of a. 

If m = 0, so that o = {v} is a vertex, ^'{v) = (p{v) and we put w^^y = v. Let 
m > and suppose p^ is constructed on simplices of lower dimensions. Let r be 
a minimal face of a such that /i(r) = Consider two cases. 

(i) If T ^ a, then Wr and ip^{wr) are defined by the induction step. We put 
Wa = Wr- Since a is convex, any x in the interior of cr is on a line segment 
joining Wa to a uniquely defined y{x) on the boundary of a. Since ip^{y{x)) 
is defined by the induction step, we extend p'' to [wa, y{x)] by linearity. 
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Figure 2: The linear (dashed line) and axis-wise linear (continuous line) interpo- 
lations of a function ip defined on the vertices of the simplex a = [vo,Vi] with 
values in R^. 

(ii) If r = a, then let Wa be the barycenter of a and put (p^{wa) = /^(o"). Again, 
any x ^ w^'m the interior of a is on a line segment joining to a uniquely 
defined y(x) on the boundary of a and we proceed as before. 

The property (a) follows from the fact that /i(r) ^ //(a) when t < a, and from 
the linearity on joining segments. The property (b) is clear from the construction. 
By routine arguments from convex analysis, the point y{x) on the boundary of a is 
a continuous function of x G cr \ {wa}, and the constructed function is continuous 
on a. Since we proceeded by induction on the dimension of a, the definitions on 
any two simplices coincide on their common face, so t/?^ extends continuously to 
K. The property (b) implies that if A; = 1, and in certain cases of vector valued 
functions, if'' is equal to Tp, namely: 

(c) if" is piecewise linear on each simplex a. In addition, if Wr is a vertex of r 
for each t < a, then it is linear on a. 

The difference between the piecewise linear and the axis-wise linear interpola- 
tions Tp and (f^ is illustrated in Figure[2lfor a 1-simplex a = [vq, vi] and a function 
(/? defined on vertices. 

Lemma 3.2 The following statements hold: 

(i) For any a G M^ Ka C K^-^^^- 

(ii) Let a E IC and a G Mf'. IfcrH K^-^^^ ^ 0, then a has at least one vertex in 

Proof: {i) Let a G /Cq. It is clear from ^ that /i(cr) -< a. It follows from the 
property (a) in the definition of (f" that a C K^-'^a- 

(ii) We follow the induction steps in the construction of ip^. If dim{a) = 0, a 
is a vertex and there is nothing to prove. Let dim(cr) = m > and suppose the 
statement is proved for lower dimensions. Let x G cr fl K^-^^^- If a; = Wa, then 
Wa G K^-^^a- By the property (a) of cp'', all a is in Ka- If x 7^ Wa, then x is on a 
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line segment joining the point w„ of a with a point y{x) of an (m — l)-simplex 
T < a, where is defined by the induction hypothesis. We know that cp^ 

is extended linearly to the line segment [wa,y{x)]. Also, (p^{y{x)) ^ (/^^(wcr) by 
the property (a). Hence (p^{y(x)) ^ f^^ix) ^ a. It remains to use the induction 
hypothesis for y{x) and r to deduce that r has a vertex in Ka- □ 

Theorem 3.3 For anj a G M^, JiTq, ?'5 a strong deformation retract of K^^^^^. 
Consequently, the inclusion Ka ^ K^^^^ induces an isomorphism in homology. 

Proof: Note that K^-^^^ is contained in a union of simplices a E JC such that 

cr^-^a := a n K^--<a ^ 0. (3) 

Given any such a, consider the simplex cTq, defined as the convex hull of the set of 
vertices v of a such that ip(v) ^ a. By the hypothesis on a and by Lemma IX2F »"). 
cTq 0. Given any cr G /C for which a^-^^^ 0' we shall define a strong 
deformation retraction 

Ha ■ Cr^^^a X [0, 1] -> (T^-^Q, 

with r = 1) being a retraction of (T<^^^a onto cTq,. 

The construction goes by induction on the dimension m of cr following the 
induction steps in the construction of the function . If dim(cr) = 0, cr is a vertex 
and there is nothing to prove. Now let m > 0. Suppose that the deformation 
retraction Hj. : t^^^q, x [0, 1] — )■ r^^^^ is defined for simplices r of dimension 
m' < m with t^^^q, 7^ in such a way that HT-{x,t) = x for any {x,t) G 
r^T^Q, X [0, 1/2™], and the values of y?^ on Hr{x,t) are decreasing with t. By 
"decreasing" we mean the weak inequality This hypothesis guarantees that 
the deformation has values in the set r^^^Q,. 

Let X G cr^^^Q. If a: is on a boundary of a, we define Ha{x,t) = Hr{x,t), 
where r is the smallest face of a containing x and Hr is defined by the induction 
hypothesis. Suppose x is in the interior of cr. Let Wa- and y{x) be the points 
identified in the definition of (/? . Note that, if (/? (wq.) ^ en, then cr^^-<a = a, 
hence the deformation must be defined as the identity map for each t: 

Ha{x,t) := X for all {x,t) e a x [0,1]. 

Therefore, we may suppose that Wa- ^ cr^^^o- Consider the smallest face r of cr 
containing y{x). Since y{x) is on the boundary of cr, r is a proper face of a of 
dimension, say, m' < m. By the construction of (/? , 

V^{.y{x)) < ip{x) < ip'iwa). (4) 

Since ^p'{x) ^ a, we get (p^{y{x)) ^ a so y{x) G t^^^q 7^ 0. By the induction 
hypothesis, a deformation retraction 

Hr ■ r^-'^a X [0, 1] T^-'-<a C CT^-^q 
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is defined so that the values of ^p" on Ht-{x, t) decrease with t, and Ht-{x, t) = x 
fort G [0,1/2™']. 

For any t G [0, 1] and for x in the interior of a we define 



It is easily checked that H„{x, 1/2™) = x, 1/2™"^) = y{x). Since ip^ 

is linear on [w^, vix)], the inequality © implies that the values of y?^ on H^{x, t) 
decrease with t. 

Thus we have defined both when x is on the boundary of a and when it 
is in the interior of a. By construction, for every x G cr^^^^a, H^(x,0) = x, 
and Hcr{x, 1) belongs to ctq, and moreover, for every x G cTq,, H^{x, 1) = x. In 
order to conclude that Ha^ is a deformation retraction of a^-^^a onto a a we must 
prove that is continuous. The continuity at a given point (xq, to) with xq in the 
interior of a follows from the continuity of y{x) in x. The continuity at (xq, to) 
with xq on the boundary of a follows from the condition that Hr{x, t) = x for any 
t G [0, 1/2™ ] and from the induction hypothesis. 

In order to continuously extend H„ to a deformation 



it is enough to prove that, given two simplices ai and a2 intersecting K^-^^^ and 
r = cTi n (T2, the maps and agree at any x G t^^^^q. It is clear from the 
definition that H„-^{x, t) = H„^{x, t) = Hr{x, t) for x G r and for all t, provided 
that Hr is defined. But this is true, because x G t^^^^q, so this is a nonempty set. 

□ 

In the next section, we use Theorem 13.31 to show that any distance between 
rank invariants of continuous functions that has property (S) can be approximated 
by the distance between rank invariants of discrete functions. We end this section 
with another application of Theorem 13.31 of interest in itself: a theorem on the 
structure of the set of critical values of the axis-wise interpolation </? . The follow- 
ing definition generalizes the notion of homological critical value given in |[T3l 
to vector functions. In plain words we call homological critical any value a for 
which any sufficiently small neighborhood contains two values whose sublevel 
sets are included one into the other but cannot be retracted one onto the other. 
Neighborhoods are taken with respect to the norm ||a|| = maxj=i,2,...,fc \oij \ in MJ'. 

Definition 3.4 Let : A' — )■ M'^ be a continuous vector function. A value a G M*^ 
is a homological critical value of if there exists an integer q such that, for all 
sufficiently small real values e > 0, two values a', a" G M.'' can be found with 

a' ^ a ^ a", \\a' — a\\ < e, \\a" — a\\ < e, such that the map 



induced by the inclusion K^^a' ^ K^^a" is not an isomorphism. If this condition 
fails, a is called a homological regular value. 




H : K^^^^ X [0, 1] ^ K^-^ 
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Also note that, by the long exact sequence for the relative homology (see e.g. 
ETl Chapter 9]) the critical value definition is equivalent to the following condi- 
tion on the graded relative homology 

For any j = 1, 2, . . . , /c and a vertex v G V(/C), consider the hyperplane of M*^ 
given by the equation aj = </3 ) and a positive closed cone Cj{v) contained in 
it, given by the formula 

Cj{v) := {a G M'' I aj = j{v) and > (f" i{v) for alH = 1, 2, . . . , k). 



Theorem 3.5 The set ofhomological critical values ofip is contained in the finite 
union of the described cones, namely, in the set 

C := \]{Cj{v) I V G V(/C) and] = 1, 2, . . . , A;}. 

Proof: Consider any a ^ C. We need to show that a is a homological regular 
value. Since C is a closed set, an e > exists such that the set Q{a,e) = {(3 G 
\\\a- I3\\< e} does not meet C. If ||a - (3\\ < e, then 

Kfs = K^. (5) 

Indeed, if this were not true, the segment joining a and (3 should contain a point 
of C, against the choice of e. 

Now, let us assume that a' < Lp{v) < a" , \\a — a'\\ < e and \\a — a"\\ < 
e. It follows from equation ([5]) and from Theorem 13.31 that the inclusions i' : 
Ka = Ka' K^-^^^i and i" : = K^" K^-^^^u induce isomorphisms in 
homology. The inclusion : K^-^^^i K^-^^^,, can be written as = 

i" o r', where r' is the retraction homotopically inverse to i' . By the functoriality 
of homology, H^{i^"''°'"'>) = H^,{i") o H^^ijr'), hence it is also an isomorphism, g 

For the sake of visualization, in Figure [3] the set C is shown in a simple case. 
From the formula for C, we instantly get an analogy of a well-known result 
from differential geometry [|24l . 

Corollary 3.6 The set ofhomological critical values ofip^ is a nowhere dense set 
in R^. Moreover its k-dimensional Lebesgue measure is zero. 

4 Approximation of distances between rank invari- 
ants 

The goal of this section is to show that shape comparison by persistent homology 
of vector functions is numerically stable. In this passage from real (continuous) 
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Figure 3: The set C defined in Theorem 13.51 is the union of closed cones with 
vertices at the values taken by tp. The set A introduced in Proposition 14.61 is a 
finite set whose elements are the highlighted points. 



objects to their discretizations, the approximation error does not grow to be much 
larger when we compute the distance D between the rank invariants of </? and '0 
instead of the distance D between the rank invariants of / and g (Theorem 14.51) . 
To this end, the stability property (S) of D defined in Section lXTl in the continuous 
setting is crucial. 

A description of this approximation procedure in concrete examples together 
with experiments exploiting the numerical stability of the comparison by persis- 
tent homology will be given in Section [51 

We end the section by showing that the set of homological critical values, al- 
though uncountable, admits a finite representative set. 

We start from the following approximation lemma. It may happen that D (p^, p^) 
and D (p^^, p^^) are equal to +oo. In the case of the matching distance, this occurs 
when H^{K) ^ H^{L). In such a case, we adopt the convention oo — oo = 0. 

Lemma 4.1 Let (p : K ^ M'^, il) : L ^ Ml^ be two continuous measuring func- 
tions on the carrier of complexes K and C. For any e > 0, there exists 5 > such 
that if 

max{diam(T | cr G /C or cr G £} < 5 (6) 

then 

|D(p^,p^) -D(p^^,p^^)| < e. (7) 

Proof: Since K and L are compact, if), Lp^ , and ip" are uniformly continuous. 
Hence for any e > 0, there exists 5 > such that if ^ is satisfied then 

max{diam(^((T) | a G /C} < e/4 (8) 

and the same inequality holds for ip" , and ilT . The diameters of (^(a) and a are 
measured with respect to the maximum norm in the respective ambient spaces. 
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Since (/? is the restriction of </? to the vertices, and (/? interpolates </? on the vertices, 
given any x G cr G /C, and any vertex u of cr, from ([8]) we get 

ll^(x) - < ll^(x) - v{v)\\ + Mv) - ^\x)\\ < e/2. (9) 

Hence, by the choice of the maximum norm in M*^, ||(^ — ||oo < e/2. By the 
same arguments, H^/) — ^/^ ||oo < e/2. By the stability property (S) of D , 

D(p^,p^) < D(p^,p^^) + D(p^^,p^,^) +D(p^-,p^J 

< ll'^-V'loo + D(p^^,p^,^) + \\lp^ -ljj\\oo 

< D (p^-,p^-) + e. 

Reversing the roles of (f, ip and ^p^, ilT, we get D (p^^ p^^) < D (p<^, p^) + e and 
the conclusion follows. □ 

Knowing Lemma |4~T1 we now turn our attention to computing D (p^^, p^^). 

The following definition sets the notation for the rank invariant of the simpli- 
cial complex filtration obtained from a discrete map Lp. Next, we show that this 
definition gives a rank invariant coinciding with the rank invariant of the contin- 
uous function ip" . Thus it is a first step in the passage from the stability of rank 
invariants for continuous functions to that of discrete ones. Moreover, this defini- 
tion is the one which we use to implement the reduction algorithm of |l5l in our 
computations in Section [51 

Definition 4.2 Consider the discrete map : V(/C) — j- M.^ defined on vertices of 
a simplicial complex /C. The q'th real space variable rank invariant or, shortly, 
q'th real rank invariant of tp is the function p^ : — N defined on each pair 
(a, (5) G as the rank of the map 

: H,{K^) ^ H,{Kp) 
induced by the inclusion map : Kp on simplicial sublevel sets. 

Theorem 4.3 Given any discrete function ip : V(/C) — )■ on the set of vertices 
of a simplicial complex JC and its axis-wise interpolation p)" , we have the equality 
of q'th real rank invariants 

Proof: Consider any (a, /3) G A';, the inclusion maps z^"'^) : K^-^^^ ^ -^</3"</3' 
and : K^^ )■ K^. Theorem 13.31 implies that for every q E 'L'^Q have the 

following commutative diagram 
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where the vertical arrows are the isomorphisms induced by the corresponding 
retractions. Thus ranking (i^"''^)) = rankiJq(j*^"'^)). 

□ 

In the sequel, we will once again use p^p to refer to real rank invariants of arbi- 
trary order. In conclusion we obtain that the distance between the rank invariants 
of two measuring (or interpolation) functions can be approximated using only the 
corresponding simplicial sublevel sets. 

Corollary 4.4 Let : K ^ M'^, ^Ij : L ^ Ml' be two continuous measuring 
functions on the carriers of complexes /C and C and let ip : V(/C) — ?■ M*', ip : 
V(£) — )■ M*^ be the discretizations of(p and ip on the sets of vertices ofK, and C, 
respectively. For any e > there exists 6 > such that if 

max{diamcr \ aE}CoraEC}<5 

then 

|D(p^,P^) -D(p^,p^)| < e. (10) 

Proof: Immediate from Lemma ITTI and Theorem 14. 3 [ □ 
We are now ready to give the main result of this section. 

Theorem 4.5 Let X and Y be homeomorphic triangulable topological spaces, 
and let f : X ^ M*^, g : Y be continuous functions. Let {K, (p) and {L, ijj), 

with K and L carriers of complexes /C' and CJ , and (p : K M'^, ^ : L — )■ R*^ 
continuous measuring functions, approximate {X, f) and {Y, g), respectively, in 
the following sense: For a fixed e > 0, there exist a homeomorphism ^ : K ^ X 
with ||<^ — / o^lloo < e/4 and a homeomorphism ( : L ^ Y with W'lp — g oC\\qo < 
e/4. Then, for any sufficiently fine subdivision K ofK' and C of C, 

|D(p/,Pg) -D(p^,p^)| < e, 

(p : V(/C) — )■ M'^, Tp : V(£) — 7- M'^ being restrictions of (p and ip on the set of 
vertices ofK, and C, respectively. 

Proof: By the triangle inequality 

D (p/, pg) < D (p/, p/oc) + D (p/oc, p^) + D (p^, p^) + D (p^, pgoc) + D {pgoc, Pg). 

Since p/ = p/o? and pg = p^oc, we have D (p/, pj^^) = and D (p^oc, Pg) = 0. 
Moreover, by the stability property (S), since ||<^ — / o ■Clloo < e/4 and — g o 
CIloo < e/4, we have D (p/o^, p^) < e/4 and D (p^, pgoc) < e/4. Therefore, 

D(p/,p,) < D(p^,p^) + e/2. 

By Corollary 14.41 there exists 6 > such that, if JC and £ are subdivisions of JC' 
and C with maxjdiam a \ a E IC or a E C} < 5, then D (p^, p^) < D (p<^, p^) + 
e/2. In conclusion we have proved that D (p/, pg) < D (p^, p^,) + e. 
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Reversing the roles of /, g and ip, ip, we get D {p^, p^) < D {pf,pg) + e, 
yielding the claim. □ 

We turn now to the question of the structure of the critical set of (/? . Recall 
from the previous section that when A; > 1, the set of homological critical values of 
a function on K with values in M'^ may be an uncountable set, although contained 
in a nowhere dense set C by Theorem 13 .51 However, the family {Ka}am'' of 
subcomplexes of /C is finite. Thus there exists a finite representative set A C M*^ 
for C as the following proposition states. 

Proposition 4.6 For any a E C, there exists A in 

A = {AgC|Vj = 1,2,...,A;, 3 G V(/C) ■\j = ^j{v)} 

such that Ka = K\. 

Proof: Since a E C, V(/Cq,) ^ and there exists j such that Uj = ^j{vj) 
for some vj G V(/C), and > ^i{vj), for 1 < i < k. For each i ^ j, let us 
take a vertex Vi G V(/Cq) such that (pi{vi) > (pi{v) for every v G V(a). Now 
we set A = (Ai, . . . , A^), with \j = Lpj(vj). By construction A belongs to A. 
Furthermore, it holds that Ka = Kx. Indeed, obviously, Kx C K^- Moreover, 
for every v G Ka, by definition of Vi it holds that Lpi(yi) > ^Pi{v) for 1 < i < k. 
Equivalently, Aj > '^i{v) iox \ < i < k, implying that v G Kx- □ 

The structure of set A is visualized in Figure [3l The previous proposition 
prompts for the following definition. 

Definition 4.7 Consider the discrete map ip : V(/C) — j- M'^ defined on vertices of 
/C. The discrete rank invariant of is the restriction of the real rank invariant 
to the finite domain A^ := n (A x A). 

Definition l4.7l gives a discrete rank invariant which is similar to the one defined 
in [HI 111, except for the fact that we are using a different homological structure. 

Defining a distance D directly on the basis of A;-dimensional rank invariants 
would be a task impossible to accomplish. Even when A; = 2, a pair of complexes 
/C and C with an order of thousand vertices would result in computing ranks of 
millions of maps induced by inclusions. This motivates the one-dimensional re- 
duction method described in the next section to compute the matching distance. 

5 Algorithm and experimentation 

For experimentation purposes, we now fix the distance between rank invariants 
that we will use to be the matching distance Dm defined in [[5l. 

The one-dimensional reduction method presented in [[5]| to compute the match- 
ing distance consists of applying the one-dimensional rank invariant along the 
lines t ^ h + tl parameterized by t and determined by pairs of vectors (/, h) in 
a chosen grid in M''' x M^, where h is an initial point and / directs the line. It is 
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assumed that all components of / are positive, and that / ■ 1 = 1, 6 ■ 1 = 0, where 

I = (1, 1, . . . , 1). For all (a, j3) E A^, there exists a unique such pair, which 
will be called linearly admissible pair or simply admissible pair, the set of which 
will be denoted Ladnik. Also denote by g^j;^^ : — > R and h^^^-^ : L — )> M the 

one-dimensional functions given by (yf^^-g^ (x) = ||(</?^(x) — h)/l\\ dMd h0-^{x) = 

II (^Z' (x) —h)/l\\, where the division is coordinate- wise and the norm is the supre- 
mum norm. For ease of notation, the pair (/, h) may be left out of g and h if it is 
unambiguous. By [|5l Lemma 1], if a = 6 + si, then 

This and Theorem 14. 31 implies 

Corollary 5.1 Consider {a, (5) = (b + sl,b + tl) G A'^, for some {s,t) G A^. 
Then 

The above theorem shows that it is legitimate to apply the reduction method of 
|I51 to simplicial sublevel sets. Following [5, Definition 11], we define the multi- 
dimensional matching distance between the rank invariants and p^ to be 

Dm (P^, P^) = _ sup . min^ k Dm (Pg^^-g, , Ph(,-gj)- 

In this section, the value minj=i,,,^A: k Dm {pg — iPg will be denoted dm (P(/j, piii) 
or dm(ff;)(p<^, Pi/i) ^nd referred to as the rescaled one-dimensional matching dis- 
tance. The computational problem is, given a threshold value e > 0, com- 
puting an approximate matching distance Dm (p,^, p^) on a suitable finite subset 
A C Ladnik such that 

Dm(P<^,P^) < Dm (p<^,P^) < Dm(p<^,Pv) + ^- (H) 

5.1 Algorithm 

Our algorithm's inputs consist of lists of simplices of K, and C of highest di- 
mension together with their adjacency relations and vertices, and of the values of 
normalized measuring functions Lp : V(/C) — and ^ : V(£) — t- M^, as well as 
a tolerance eQ In our computations, we have confined ourselves to the case where 
K, and C are triangular meshes. Its output is an approximate matching distance 
Dm (P(p5 Pv)- compute the one-dimensional persistent homology on admissi- 
ble pairs, we use the persistent homology software JPlex |fT6l . By default, JPlex 
computes the persistent Betti numbers over Zn of a discretely indexed filtration of 
simplicial complexes. We build this filtration by adding simplices in the following 

'Due to the finite precision of computer arithmetic, the codomain of the functions Lpi and ■4>i is 
in reality 10"^ 1 rather than R. In our computations we tended to use p ~ Q, that is, a precision 
of up to six digits after the decimal point. 
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recursive way. We first order the values attained by the one-dimensional measur- 
ing function g in increasing order, {gi, . . . , g^}- A finite filtration {/Ci, . . . , /Cat} 
is then built by inserting simplices a G /C. If a = {v}, where f is a vertex, we put 
{v} into JCi if g{v) < gi. Otherwise, a G /C/ if all its vertices are in JCi. Similarly, 
the function h is used to build a finite filtration {Ci, . . . , Cm} using simplices of 
C. 

The set of admissible pairs Ladm2 is the set of quadruples (a, 1 — a,b, —b) G 
X such that < a < 1. As described in [|2l Remark 3.2], it is possible to 
avoid computation of the one-dimensional matching distance over a large portion 
of Ladm2. Since the functions (p and ^ are normalized, C = max{\\Lp\\, Wi'W} = 
1. Let Ladml be the set of admissible pairs such that \b\ < 1. Then, to compute 
the maximal value of Dm {pg, p^j over Ladm2\Ladm2, it is sufficient to consider 
the two admissible pairs (a, 1-a, b, -b) = (1/2, 1 /2, 2, -2) and (1/2, 1/2, -2, 2). 
The details can be found in [|2l. 

It follows from a generalization of the Error Bound Theorem (|l2l Theorem 3.4] 
and ifTOll ) to persistent homology of arbitrary order that if for (/, b) and (/', b') G 
Ladm2, \\ (/, b) — (/', b') \\ < 5, then for normalized functions if and %p 

I dm(;;b)(P¥',Pv) ^ dm(;?_g/)(p<^,P^)| < 185. 

This suggests that in order to satisfy Equation (fTTI) . it suffices to choose admissible 
pairs (/, 6) G Ladrn'^ at a distance within e/9 of each other, guaranteeing that 
every member of Ladm^ is within e/18 of a tested pair. In practice, our algorithm 
is reminiscent of the grid algorithm shown in Section 3 of O, in the sense that 
we take pairs at a distance of 1/2^ of each other with N sufficiently large. We 
observe that the set Ladm2 is in bijective correspondence with {(a, fe) G | a G 
(0, 1), b G M}, and so we will speak of computing dm(Pi^, p^.) at apoint P = (a, b) 
of the preceding set. The lattice of points on which we compute this rescaled 
matching distance is chosen as follows: choose N E N such that 1/2^ < e/18, 
and choose Pij = (a^, bj),i = 0, . . . , 2^ - 1, j = 0, . . . , 2^+^ - 1 such that 
a, = (2z + l)/2^+i, b, = l- (2j + l)/2^+\ 

5.2 Examples of topological aliasing 

Our experimentations have been made on triangular meshes of compact 2D sur- 
faces. In doing so, the influence on experimental results of the concept of topo- 
logical aliasing discussed in Section [3] became apparent. Namely, we used our 
algorithm to compare in a pairwise manner 10 cat models, a selection of which 
is found in Figure IH We used for Lfi and ^i, i = 1,2, the following functions. 
Assume that the model /C is such that its vertex set V(/C) = {vi, . . . , tv} and 
compute the following principal vector: 

Er=iK-c)ii^»-cii2 

\-^n II 112 ' 

where c is the centre of mass of K defined by taking the weighted average of 
the centres of each triangle. Let d be the line passing through c having w as its 
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direction vector, and let vr be the plane passing through c having tv as its normal 
vector. We defined 

dist{vi,d) 

ipiiVi) = 1 — — — 

maxj=i_...^„ dist^Wj, a) 

and 

dist(?;i,7r) 
¥^2 = 1 jr-j r, 

where dist(t;, d) and dist(t;, it) are defined in the usual way, as the minimal Eu- 
clidean distance between v and the points on d or vr. The functions 'tpi and tp2 were 
defined similarly using the model £. We then repeated the same procedure on the 
barycentric subdivisions of the models, with the value of the function at the new 
vertices defined using the linear interpolant. We found out that in this case the 
computed matching distance did not always yield the same result as when using 
the original unsubdivided models. However, replacing the linear interpolant by 
the axis- wise linear interpolant allowed us to retrieve the same results. 

We show in Table [T] a selected subset of our results, and in Figure |4] images 
of the five models for which results are shown. The first two rows of numbers 
in each table represent the ID matching distance computed using each of (pi and 
tpi for 2 = 1,2 respectively, while the last three rows represent the approximated 
2D matching distance computed for three different tolerance levels e. The col- 
umn named "Nonsub" shows the distances computed on the original model, while 
"Linear" and "Axis-wise" show those computed on the subdivided models with re- 
spectively the linear and axis-wise linear interpolants. "Diff" and "% Diff" show 
the difference and relative difference between the matching distance results for 
the unsubdivided models and subdivided models with linear interpolation. 

We can see that while the matching distance computed using the axis-wise lin- 
ear interpolant is, for every tolerance level, equal to the matching distance between 
the original models, the matching distance computed using the linear interpolant 
can be quite different, and this both using 0th- and Ist-order persistent homology. 
However, this phenomenon is only seen when computing the 2D matching dis- 
tance: the ID matching distances (the numbers in the first two rows of each table) 
are always the same. Given that in our context topological aliasing can only be 
observed when using multi-measuring functions, this follows our expectations. 
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Figure 4: The models catO, catO-tranl-1, catO-tranl-2, cat0-tran2-l and catO- 
tran2-2 are shown along with the values of Lpi and Lp2- Models are courtesy of 
the authors of [2J. 
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catO vs. catO-tranl-1 




Nonsub 


Linear 


Axis-wise 


Dili 


% Dili 




0.031129 


0.031129 


0.031129 


0.000000 


0.000000 




0.039497 


0.039497 


0.039497 


0.000000 


0.000000 


Hi 


0.039497 


0.039497 


0.039497 


0.000000 


0.000000 




0.046150 


0.039497 


0.046150 


-0.006653 


-16.844317 




0.046150 


0.040576 


0.046150 


-0.005574 


-13.737185 




0.118165 


0.118165 


0.118165 


0.000000 


0.000000 




0.032043 


0.032043 


0.032043 


0.000000 


0.000000 


Ho 


0.194217 


0.177001 


0.194217 


-0.017216 


-9.726499 




0.224227 


0.203102 


0.224227 


-0.021125 


-10.401178 




0.225394 


0.207266 


0.225394 


-0.018128 


-8.746249 





catO-tranl-2 vs. cat0-tran2-l 




Nonsub 


Linear 


Axis-wise 


Diff 


% Diff 




0.017272 


0.017272 


0.017272 


0.000000 


0.000000 




0.026101 


0.026101 


0.026101 


0.000000 


0.000000 


Hi 


0.026101 


0.028686 


0.026101 


0.002585 


9.903835 




0.034314 


0.028686 


0.034314 


-0.005628 


-19.619327 




0.034314 


0.029188 


0.034314 


-0.005126 


-17.562012 




0.182985 


0.182985 


0.182985 


0.000000 


0.000000 




0.018951 


0.018951 


0.018951 


0.000000 


0.000000 


Ho 


0.192872 


0.188365 


0.192872 


-0.004507 


-2.392695 




0.207480 


0.202844 


0.207480 


-0.004636 


-2.285500 




0.208451 


0.204511 


0.208451 


-0.003940 


-1.926547 














cat0-tran2-l vs. cat0-tran2-2 




Nonsub 


Linear 


Axis-wise 


Diff 


% Diff 




0.022001 


0.022001 


0.022001 


0.000000 


0.000000 




0.034288 


0.034288 


0.034288 


0.000000 


0.000000 


Hi 


0.034288 


0.034288 


0.034288 


0.000000 


0.000000 




0.045545 


0.035702 


0.045545 


-0.009843 


-27.569884 




0.045545 


0.037061 


0.045545 


-0.008484 


-22.891989 




0.095677 


0.095677 


0.095677 


0.000000 


0.000000 




0.032966 


0.032966 


0.032966 


0.000000 


0.000000 


Ho 


0.178776 


0.182322 


0.178776 


0.003546 


1.983488 




0.202770 


0.196977 


0.202770 


-0.005793 


-2.940952 




0.212733 


0.208097 


0.212733 


-0.004636 


-2.227807 



Table 1: The approximated matching distance computed by our algorithm for 
three decreasing tolerance values (e = 9/8,9/16 and 9/32) is shown for a few test 
cases (unsubdivided, subdivided with ^ and subdivided with ), with 0th- and 
Ist-order rank invariants. 
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